!>\brief
!! This module collects the computation of the local matrices for the
!! LDG-H method.
!!
!! \n
!!
!! Essentially, we compute here the element matrices and right hand
!! sides, before any static condensation. These terms include Neumann
!! type boundary conditions, while ignore Dirichlet type boundary
!! conditions (which enforcement in essential form must be done
!! elsewhere).
!! \note There are no consistency checks on the coefficients returned
!! by the \c coeff* functions of \c mod_testcases, which thus have to
!! be correct. In particular, the diffusion tensor must be symmetric.
!!
!! \par Internal barriers
!!
!! If \c coeff_alpha is associated, the local problems are treated as
!! Robin problems, with a finite jump between the element primal
!! variable and the hybrid variable. This requires the following
!! changes:
!! <ul>
!!  <li> The Lagrange multiplier is introduced in the weak form of the
!!  constitutive equation with the formal substitution
!!  \f{displaymath}{
!!    u|_{\partial K} \mapsto \lambda + \alpha\left(\hat{q}\cdot
!!    n\right).
!!  \f}
!!  <li> The numerical flux is defined as
!!  \f{displaymath}{
!!    \hat{q} = \frac{1}{1+\alpha\tau} q  +
!!    \frac{\tau}{1+\alpha\tau}(u-\lambda)n.
!!  \f}
!!  Notice that this choice is consistent, since
!!  \f{displaymath}{
!!    \hat{q} = q \Longleftrightarrow -\alpha\left(q\cdot n\right)
!!    + u = \lambda.
!!  \f}
!!  Notice as well that \f$\hat{q}\f$ is (formally) well defined for
!!  \f$\alpha\in[0\,,\infty]\f$.
!! </ul>
!! \note In presence of a fracture, fhe flux balance low should
!! include also the along-fracture flow; more precisely, the weak form
!! of the Darcy problem along the fracture yields the generalized
!! conservativity condition for neighboring elements. This is not
!! taken into account in this module, however, since the local
!! problems define only the outgoing flux \f$\hat{q}\f$. Enforcing the
!! conservativity condition is something that is performed during the
!! assembly of the global system, and that is when the along-fracture
!! flow must be included.
!! \note The case \f$\alpha=0\f$ coincides with the usual one, where
!! the hybrid variable provides a Dirichlet datum for the local
!! problem.
!<----------------------------------------------------------------------
module mod_ldgh_locmat

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_linal, only: &
   mod_linal_initialized, &
   invmat_chol
 
 use mod_master_el, only: &
   mod_master_el_initialized, &
   t_me

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_grid, only: &
   mod_grid_initialized, &
   t_e, affmap, t_ddc_grid

 use mod_bcs, only: &
   mod_bcs_initialized, &
   t_bcs,                     &
   t_b_v,   t_b_s,   t_b_e,   &
   p_t_b_v, p_t_b_s, p_t_b_e, &
   b_dir,   b_neu,   b_ddc
 
 use mod_tau, only: &
   mod_tau_initialized, &
   t_eltau

 use mod_testcases, only: &
   mod_testcases_initialized, &
   coeff_diff,  &
   coeff_adv,   &
   coeff_re,    &
   coeff_f,     &
   coeff_neu,   &
   coeff_rob,   &
   coeff_alpha

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_ldgh_locmat_constructor, &
   mod_ldgh_locmat_destructor,  &
   mod_ldgh_locmat_initialized, &
   ldgh_locmat, &
   t_bcs_error

 private

!-----------------------------------------------------------------------

! Module types and parameters
 
 !> This type is used to indicate that an error has occurred in the
 !! computation of the boundary condition
 type t_bcs_error
   logical :: lerr = .false.
   character(len=1000) :: message
 end type t_bcs_error

! Module variables

 ! public members
 logical, protected ::               &
   mod_ldgh_locmat_initialized = .false.
 ! private members
 real(wp), allocatable :: ddkr(:,:), mmkr(:,:)
 character(len=*), parameter :: &
   this_mod_name = 'mod_ldgh_locmat'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_ldgh_locmat_constructor(base)
  type(t_base), intent(in) :: base

  integer :: l, i, j, is
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_linal_initialized.eqv..false.)    .or. &
       (mod_master_el_initialized.eqv..false.).or. &
       (mod_base_initialized.eqv..false.)     .or. &
       (mod_grid_initialized.eqv..false.)     .or. &
       (mod_bcs_initialized.eqv..false.)      .or. &
       (mod_tau_initialized.eqv..false.)      .or. &
       (mod_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_ldgh_locmat_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   allocate( ddkr(base%pk,base%mk) )
   ddkr = 0.0_wp
   do l=1,base%m
     do i=1,base%pk
       do j=1,base%mk
         ddkr(i,j) = ddkr(i,j) - base%wg(l) *         &
           dot_product(base%o(:,j,l),base%gradp(:,i,l))
       enddo
     enddo
   enddo
   !  boundary contributions
   do is=1,base%me%d+1
     do l=1,base%ms
       do i=1,base%pk
         do j=1,base%mk
           ddkr(i,j) = ddkr(i,j) + base%wgb(l,is) * &
             base%no(j,l,is)*base%pb(i,l,is)
         enddo
       enddo
     enddo
   enddo

   allocate(mmkr(1,base%pk))
   mmkr = 0.0_wp
   do l=1,base%m
     do j=1,base%pk
       mmkr(1,j) = mmkr(1,j) + base%wg(l) * base%p(j,l)
     enddo
   enddo

   mod_ldgh_locmat_initialized = .true.
 end subroutine mod_ldgh_locmat_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_ldgh_locmat_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_ldgh_locmat_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate(ddkr,mmkr)

   mod_ldgh_locmat_initialized = .false.
 end subroutine mod_ldgh_locmat_destructor

!-----------------------------------------------------------------------
 
 !> Compute the local matrices and right hand side.
 !!
 !! The local matrices correspond to the following terms:
 !! \f{displaymath}{
 !!  A^K_{ij} = \int_K \omega_j\cdot \mu^{-1} \omega_i\,dx +
 !!  \int_{\partial K} c_2 \omega_{n_j} \omega_{n_i} d\sigma, \qquad
 !!  B^K_{ij} = -D^K_{ji} - \int_K \varphi_j\mu^{-1}a\cdot\omega_i\,dx,
 !!   \qquad
 !!  C^K_{ij} = F^K_{ji},
 !! \f}
 !! \f{displaymath}{
 !!  D^K_{ij} = \int_K \nabla\cdot\omega_j\,\varphi_i \,dx -
 !!  \int_{\partial K} c_4 \omega_{n_j}\varphi_i\, d\sigma,
 !!   \qquad
 !!  E^K_{ij} = \int_K \sigma\varphi_j\varphi_i\,dx
 !!             +\int_{\partial K}c_3\varphi_j\varphi_i\,d\sigma,
 !!   \qquad
 !!  H^K_{ij} = -G^K_{ji},
 !!   \qquad
 !!  f^{K,2}_{i} = \int_K f\varphi_i\,dx,
 !! \f}
 !! \f{displaymath}{
 !!  F^K_{ij} = \int_{\partial K}c_1\omega_{n_j}\eta_i\,d\sigma,
 !!   \qquad
 !!  G^K_{ij} = \int_{\partial K}c_3 \varphi_j\eta_i\,d\sigma,
 !!   \qquad
 !!  L^K_{ij} = -\int_{\partial K}c_3\eta_j\eta_i\,d\sigma \,+\, bcs,
 !!   \qquad
 !!  f^{K,3}_{i} = bcs.
 !! \f}
 !! The coefficients \f$c_*\f$ are defined as
 !! \f{displaymath}{
 !!   c_1(\alpha,\tau) = \frac{1}{1+\alpha\tau}
 !! \f}
 !! \f{displaymath}{
 !!   c_2(\alpha,\tau) = \alpha c_1(\alpha,\tau) =
 !!   \frac{\alpha}{1+\alpha\tau}
 !! \f}
 !! \f{displaymath}{
 !!   c_3(\alpha,\tau) = \tau c_1(\alpha,\tau) =
 !!   \frac{\tau}{1+\alpha\tau}
 !! \f}
 !! \f{displaymath}{
 !!   c_4(\alpha,\tau) =
 !!  \left\{
 !!   \begin{array}{c}
 !!    \tau c_2(\alpha,\tau) \\\
 !!    \alpha c_3(\alpha,\tau)
 !!   \end{array}
 !!  \right\}
 !!  = \frac{\alpha\tau}{1+\alpha\tau}.
 !! \f}
 !!
 !! \note The optional arguments \c ddc_grid and \c me_g2l must be
 !! passed when the matrices are required in terms of the global
 !! degrees of freedom, typically during the construction of the
 !! global linear system. Otherwise, the matrices are computed in
 !! terms of the degrees of freedom of the present subdomain: this is
 !! typically what is used during the post-processing.
 !!
 !! When using a Lagrange multiplier to enforce zero mean of the
 !! primal unknown, the following matrix is also computed:
 !! \f{displaymath}{
 !!  M_{1j}^K = \int_K \varphi_j\,dx.
 !! \f}
 pure subroutine ldgh_locmat(aak , bbk , cck ,      &
                             ddk , eek , hhk , f2k, &
                             ffk , ggk , llk , f3k, &
                             base,e,tau,be,bcs_err, &
                             ddc_grid, me_g2l,      &
                             zero_mean , mmk)
  real(wp), intent(out) :: &
    aak(:,:) , bbk(:,:) , cck(:,:) ,         &
    ddk(:,:) , eek(:,:) , hhk(:,:) , f2k(:), &
    ffk(:,:) , ggk(:,:) , llk(:,:) , f3k(:)
  type(t_base), intent(in) :: base
  type(t_e), intent(in) :: e
  type(t_eltau), intent(in) :: tau
  type(p_t_b_e), intent(in) :: be
  type(t_ddc_grid), intent(in), optional :: ddc_grid
  real(wp),         intent(in), optional :: me_g2l(:,:,:)
  logical,          intent(in), optional :: zero_mean
  real(wp),        intent(out), optional :: mmk(:,:)
  type(t_bcs_error), intent(out) :: bcs_err

  integer :: l, i, j, is, shift_is, jj, ii
  real(wp) :: &
    mu(e%m,e%m,base%m), mui(e%m,e%m,base%m), alpha(base%ms),      &
    c1(base%ms,e%d+1), c2(base%ms,e%d+1), c3(base%ms,e%d+1),      &
    c4(base%ms,e%d+1), a(e%m,base%m), sigma(base%m), f(base%m),   &
    mui_b(e%m,e%d), bt_mui_b(e%d,e%d), a_mui_b(e%d), mu_o(e%d),   &
    a_mu_o, swg, swg1, swg2, xg(e%m,base%m), wgsa(base%ms,e%d+1), &
    esa(base%nk,base%ms,e%d+1)
  character(len=*), parameter :: &
    this_sub_name = 'ldgh_locmat'

   ! Precompute side weights and basis functions taking into account
   ! the proper order of the quadrature nodes (in fact, the reordering
   ! of the quadrature weights is not necessary, because the
   ! quadrature formulas are symmetric and the weights are equal for
   ! corresponding nodes, but we reorder the weights as well for
   ! clarity and consistency with the reordering of the basis
   ! functions).
   do is=1,base%me%d+1
     wgsa(:,is) = (e%s(is)%p%a/base%me%voldm1) * &
                   base%wgs(base%stab(e%ip(is),:))
     do j=1,base%nk
       esa(j,:,is) = base%e(j,base%stab(e%ip(is),:))
     enddo
   enddo
 
   ! Gaussian nodes in physical space
   xg = affmap(e,base%xig)

   ! problem coefficients
   mu    = coeff_diff(xg)
   do l=1,base%m
     call invmat_chol( mu(:,:,l) , mui(:,:,l) )
   enddo
   a     = coeff_adv(xg)
   sigma = coeff_re(xg)
   f     = coeff_f(xg)

   ! now the coefficient \alpha, defined on the element sides
   if(associated(coeff_alpha)) then
     do is=1,base%me%d+1
       alpha = coeff_alpha( affmap(e,base%xigb(:,:,is)) )
       c1(:,is) = 1.0_wp/(1.0_wp + alpha*tau%tau(:,is))
       c2(:,is) =         alpha*c1(:,is)
       c3(:,is) = tau%tau(:,is)*c1(:,is)
       c4(:,is) =         alpha*c3(:,is)
     enddo
   else
     c1 = 1.0_wp
     c3 = tau%tau
   endif

   !--------------------------------------------------------------------
   ! First equation: compute aak and the transport contribution to bbk

   ! first the volume terms
   aak = 0.0_wp
   bbk = 0.0_wp
   do l=1,base%m
     ! scaled quad weight
     swg1 = base%wg(l)/e%det_b
     swg2 = base%wg(l)
     ! mui_b = mui * b
     mui_b = matmul( mui(:,:,l) , e%b )
     ! bt_mui_b = b^T * mui * b
     bt_mui_b = matmul( transpose(e%b) , mui_b )
     ! a_mui_b = a^T * mui * b
     a_mui_b = matmul( a(:,l) , mui_b )
     do i=1,base%mk
       ! mu_o = b^T * mui * b * o
       mu_o = matmul( bt_mui_b , base%o(:,i,l) )
       do j=1,base%mk
         aak(i,j) = aak(i,j) +                     &
           swg1 * dot_product( base%o(:,j,l) , mu_o )
       enddo
       ! a_mu_o = a^T * mui * bk * o   (including swg2)
       a_mu_o = swg2*dot_product(a_mui_b,base%o(:,i,l))
       do j=1,base%pk
         bbk(i,j) = bbk(i,j) - base%p(j,l)*a_mu_o
       enddo
     enddo
   enddo
   ! now the boundary terms
   if(associated(coeff_alpha)) then
     do is=1,base%me%d+1
       do l=1,base%ms
         swg = (base%me%a(is)/e%s(is)%p%a)*base%wgb(l,is) * c2(l,is)
         do i=1,base%mk
           do j=1,base%mk
             aak(i,j) = aak(i,j) +                   &
               swg * base%no(j,l,is) * base%no(i,l,is)
           enddo
         enddo
       enddo
     enddo
   endif
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Second equation: compute ddk, eek and f2k

   ! first the volume terms
   ddk = ddkr
   eek = 0.0_wp
   f2k = 0.0_wp
   do l=1,base%m
     swg1 = e%det_b*base%wg(l) * sigma(l)
     swg2 = e%det_b*base%wg(l) * f(l)
     do i=1,base%pk
       ! ddk volume contributions are precomputed
       do j=1,base%pk
         eek(i,j) = eek(i,j) +          &
           swg1 * base%p(j,l)*base%p(i,l)
       enddo
       f2k(i) = f2k(i) +  &
         swg2 * base%p(i,l)
     enddo
   enddo
   ! now the boundary terms
   if(associated(coeff_alpha)) then
     do is=1,base%me%d+1
       do l=1,base%ms
         swg = base%wgb(l,is) * c4(l,is)
         do i=1,base%pk
           do j=1,base%mk
             ddk(i,j) = ddk(i,j) -                   &
               swg * base%no(j,l,is) * base%pb(i,l,is)
           enddo
         enddo
       enddo
     enddo
   endif
   do is=1,base%me%d+1
     do l=1,base%ms
       swg = wgsa(l,is) * c3(l,is)
       do i=1,base%pk
         do j=1,base%pk
           eek(i,j) = eek(i,j) +                   &
             swg * base%pb(j,l,is) * base%pb(i,l,is)
         enddo
       enddo
     enddo
   enddo
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Third equation: compute ffk, ggk, llk and f3k

   ! here there are only boundary terms
   ffk = 0.0_wp
   ggk = 0.0_wp
   llk = 0.0_wp
   f3k = 0.0_wp
   do is=1,base%me%d+1
     shift_is = (is-1)*base%nk
     do l=1,base%ms
       swg1 = base%wgb(l,is) * c1(l,is)
       swg2 =     wgsa(l,is) * c3(l,is)
       do i=1,base%nk
         ii = shift_is + i
         do j=1,base%mk
           ffk(ii,j) = ffk(ii,j) +              &
             swg1 * base%no(j,l,is) * esa(i,l,is)
         enddo
         do j=1,base%pk
           ggk(ii,j) = ggk(ii,j) +              &
             swg2 * base%pb(j,l,is) * esa(i,l,is)
         enddo
         do j=1,base%nk
           jj = shift_is + j
           llk(ii,jj) = llk(ii,jj) -        &
             swg2 * esa(j,l,is) * esa(i,l,is)
         enddo
       enddo
     enddo
   enddo

   ! add boundary fluxes if present
   bcs_err%lerr = .false.
   if(associated(be%p)) then
     call ldgh_boundary_flux(ffk,ggk,llk,f3k,base,be%p, &
                             ddc_grid,me_g2l, bcs_err)
   endif
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Fourth equation: enforce zero mean
   if(present(zero_mean)) then
     if(zero_mean) mmk = e%det_b * mmkr
   endif
   !--------------------------------------------------------------------

   !--------------------------------------------------------------------
   ! Now we can exploit the symmetry of the system
   bbk = bbk - transpose(ddk)
   cck =  transpose(ffk)
   hhk = -transpose(ggk)
   !--------------------------------------------------------------------

 end subroutine ldgh_locmat
 
!-----------------------------------------------------------------------
 
 !> Add the boundary fluxes.
 !!
 !! If we call this function, it means that the element has at least
 !! one boundary side. Indeed, the number of boundary sides is
 !! <tt>be\%nbs</tt>. In this subroutine we loop on the boundary sides
 !! and add the corresponding terms.
 pure subroutine ldgh_boundary_flux(ffk,ggk,llk,f3k,base,be, &
                                    ddc_grid,me_g2l, err)
  real(wp), intent(inout) :: ffk(:,:), ggk(:,:), llk(:,:), f3k(:)
  type(t_base), intent(in) :: base
  type(t_b_e), intent(in) :: be
  type(t_ddc_grid), intent(in), optional :: ddc_grid
  real(wp), intent(in), optional :: me_g2l(:,:,:)
  type(t_bcs_error), intent(out) :: err
 
  integer :: isb, isl, breg, shift_is, l, i, ii, j, jj
  real(wp) :: xgs(be%e%m,base%ms), a(be%e%m,base%ms),      &
    gamma(base%ms), hb(base%ms), an(base%ms), ap(base%ms), &
    wgsa(base%ms), esa(base%nk,base%ms),                   &
    mm(base%nk,base%nk), mt(base%nk,base%nk)

   err%lerr = .false.
   do isb=1,be%nbs
     btype: select case(be%bs(isb)%p%bc)
      case(b_dir)
       ! nothing to do in this subroutine
      case(b_neu)

       isl = be%bs(isb)%p%s%isl(1) ! local side index on element e
       shift_is = (isl-1)*base%nk
       ! Gaussian nodes in physical space
       xgs = affmap(be%e,base%xigb(:,:,isl))
       ! side weigths reordered and scaled by the side area
       wgsa = (be%e%s(isl)%p%a/base%me%voldm1) * &
                   base%wgs(base%stab(be%e%ip(isl),:))
       do j=1,base%nk
         esa(j,:) = base%e(j,base%stab(be%e%ip(isl),:))
       enddo

       ! problem coefficients
       breg = -be%bs(isb)%p%s%ie(2) ! boundary region
       gamma = coeff_rob(xgs,breg)
       hb    = coeff_neu(xgs,breg)

       btypespec: select case(be%bs(isb)%p%btype)

        case(1) ! Neumann, Hughes flux, hybrid variable

         ! advection coefficient
         a  = coeff_adv(xgs)
         an = matmul(transpose(a),be%e%n(:,isl))
         ap = 0.5_wp*(an + abs(an))
         
         ! boundary contributions to llk and f3k
         do l=1,base%ms
           do i=1,base%nk
             ii = shift_is + i
             do j=1,base%nk
               jj = shift_is + j
               llk(ii,jj) = llk(ii,jj) +                       &
                   (-ap(l)+gamma(l))*wgsa(l) * esa(j,l)*esa(i,l)
             enddo
             f3k(ii) = f3k(ii) +        &
                 hb(l)*wgsa(l) * esa(i,l)
           enddo
         enddo

        case default
         err%lerr = .true.
         write(err%message,'(A,I5,A,I5,A,I5,A)') &
           'Unknown boundary condition type "',be%bs(isb)%p%btype, &
           '" on side ',be%bs(isb)%p%s%i,', on element ',be%e%i,'.'
       end select btypespec

      case(b_ddc)
       ! The local matrices involving the side basis functions must be
       ! transformed as described in mod_ldgh_linpb.
       if(present(ddc_grid)) then
         ii = ddc_grid%gs( be%bs(isb)%p%s%i )%ni 
         if(ddc_grid%ns(ii)%id.lt.ddc_grid%id) then
           isl = be%bs(isb)%p%s%isl(1) ! local side index on element e
           shift_is = (isl-1)*base%nk

           ! get the index of the required permutation
           i = ddc_grid%ns(ii)%p_s2s
           mm = me_g2l(:,:,i)
           mt = transpose(mm)

           ! multiply the matrices
           i  = shift_is+1
           ii = shift_is + base%nk
           ffk(i:ii, : )  = matmul( mt , ffk(i:ii, : ) )
           ggk(i:ii, : )  = matmul( mt , ggk(i:ii, : ) )
           llk(i:ii,i:ii) = matmul( mt , matmul( llk(i:ii,i:ii) , mm ) )
           ! f3k is zero for internal sides
           ! nothing to do for cck: exploit cck = ffk^T, same for hhk
         endif
       endif

      case default
       err%lerr = .true.
       write(err%message,'(A,I5,A,I5,A,I5,A)') &
         'Unknown boundary condition "',be%bs(isb)%p%bc,        &
         '" on side ',be%bs(isb)%p%s%i,', on element ',be%e%i,'.'
     end select btype
   enddo
 
 end subroutine ldgh_boundary_flux
 
!-----------------------------------------------------------------------

end module mod_ldgh_locmat

